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Abstract 

We propose an efficient numerical scheme, based on the method of lines, for solving the Landau- 
de Gennes equations describing the relaxational dynamics of nematic liquid crystals. Our method is 
computationally easy to implement, balancing requirements of efficiency and accuracy. We bench- 
mark our method through the study of the following problems: the isotropic- nematic interface, 
growth of nematic droplets in the isotropic phase and the kinetics of coarsening following a quench 
into the nematic phase. Our results, obtained through solutions of the full coarse-grained equations 
of motion with no approximations, provide a stringent test of the de Gennes ansatz for the isotropic 
- nematic interface, illustrate the anisotropic character of droplets in the nucleation regime and 
validate dynamical scaling in the coarsening regime. 

PACS numbers: 64.60. Cn, 61.30.Jf, 82.20Mj, 05.70.Fh 
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I. INTRODUCTION 



Liquid crystalline phases of matter are rich in examples of subtle order parameters, com- 
plex phase behaviour, and exotic defect structures In the nematic phase of liquid crys- 
tals, broken rotational symmetry produces elastic and hydrodynamic modes and topological 
defects of integer and half- integer charge Understanding the statics and dynamics of 
nematics is difficult because the order parameter is neither a scalar or a vector, but a more 
complicated tensorial quantity constrained by symmetry and normalisation [l]]. For example, 
the simplest relaxational dynamics which follows from a time-dependent Ginzburg-Landau 
equation describing a nematic close to thermal equilibrium, entails the solution of a set of 5 
coupled non-linear parabolic partial differential equations for the independent components 
of the order parameter tensor Q jsl- For all but the simplest situations, these equations do 
not have analytical solutions and thus need to be solved numerically. The efficient numerical 
computation of the solutions of the relaxational dynamics of nematics - nematodynamics - 
is the problem addressed here. 

In this paper, we propose an efficient numerical scheme, based on the method of lines 
{4, 5], for solving the Landau-de Gennes equations for nematodynamics. The method of 
lines (MOL) is a powerful technique for discretizing initial-value partial differential equa- 
tions (PDEs). The essence of the MOL is semi-discretisation, where a PDE in spatial and 
temporal variables is discretised in the spatial variable only. This reduces the PDE to a 
system of ODEs in the temporal variable. The great advantage of the MOL is that there 
are powerful numerical methods implemented in general-purpose numerical libraries for the 
solution of systems of ODEs. The MOL also provides a great degree of freedom in discretis- 
ing space, allowing one to chose from finite-difference, finite-volume, finite-element, spectral 
collocation, or any other suitable spatial discretisation. The strategy of semi-discretisation, 
where the spatial discretisation and temporal integration are treated as separate steps, al- 
lows for considerable fiexibility, as both these operations can be optimised according to the 
demands of the problem to ensure maximum accuracy at minimum computational cost. 

Here we combine a finite-difference spatial discretisation with a Runge-Kutta temporal 
integration to solve the relaxational equations of nematodynamics. We benchmark our 
method through the study of the following three problems: (a) the isotropic - nematic 
interface and tests of the de Gennes ansatz, (b) the kinetics of the growth of single nematic 
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droplets in the isotropic phase and (c), coarsening kinetics following a quench from the 
isotropic into the nematic phase. To enable a comparison with previously published results 
our numerical results are in two dimensions, though the method and our code are applicable 
in one, two and three dimensions. 

Our results are summarised as follows. Our numerical scheme permits a stringent test of 
the de Gennes ansatz for the uniaxial nematic isotropic interface [s]. We show that biaxiality 
is absent in the interfacial region and that the hyperbolic tangent profile proposed by de 
Gennes provides an accurate fit to the numerical data. We show next that single nematic 
droplets placed in an isotropic solvent at isotropic-nematic coexistence coarsen anisotropi- 
cally j^, with initially circular droplets deforming into more ellipsoidal structures as time 
evolves. In our study of coarsening upon quenches from the isotropic phase, we observe 
the characteristic annihilation of integer and half integer charge defects (disclinations) as 
revealed by schlieren plots. For the coarsening problem, we find that order parameter cor- 
relation functions at different times can be rescaled to lie on a master curve, and that the 
coarsening exponent associated with the growing length scales is |. Our numerical results 
are in excellent agreement with analytical results where available and consistent with nu- 
merical results where such result exists. In addition, we extend previous results for each of 
the problems we study in significant ways. 

In Section [TTl we summarise the order parameter description of nematics in terms of the 
symmetric, traceless order parameter tensor Q^p- We discuss the Landau - de Gennes free 
energy functional describing the free energy associated with order parameter configurations, 
the phase diagram for the free energy and the equations of relaxational dynamics. In section 
mil we present the MOL scheme for partial differential equations and show how it may be 
applied to nematodynamics. In section IIVI we discuss each of the applications described 
above in some detail, presenting our numerical results. We conclude by describing extensions 
of the present method to other problems in the dynamics of nematic liquid crystals and 
compare our method with other methods available in the literature. 



II. RELAXATIONAL DYNAMICS 



Orientational order in t 
ric traceless tensor, Q^p [l 



le nematic phase is quantified through a second-rank, symmet- 
. The principal axes of this tensor, obtained by diagonalizing 
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Q, specify the direction of ordering. The principal values represent the strength of order- 
ing. Locally, nematic order is defined through Q(x, t) = J duf{:s., u, t) uu' = ('uu') where 
/(x, u, t) is the molecular orientational distribution function at the position x at time t 
which counts the number of molecules with long axis oriented in the direction u, X defines 
the symmetric traceless part of an arbitrary real second rank tensor X with 

X =^(X + X^)-iTr(X), (1) 

where denotes the transpose of the tensor i.e. (X"^)^^ = Xj3a and d denotes the dimen- 
sion. The average < • >= / duf{'x,u,t){-). 

The three principal axes specify the director n, the codirector 1 and the joint normal to 
these, m. Given the two principal values S and T, the order parameter can be written as 

3 11 

where a, P = x,y, z and, 

s = {p,(cose)) = {{cos'e~^-)), (3) 

T = (sin2^cos20), (4) 

are measures of alignment with — |<S'<|,0<T< 35*. 5' = | corresponds to the 
nematic phase whereas 5 = to the isotropic phase. The (polar) angle between n and u is 
6, while the (azimuthal) angle between 1 and u is 0. 

In the uniaxial nematic phase, T = and the tensor ([2]) takes the form Q = (3/2)S''nn' . 



In the biaxial phase Q = (3/2)5* nn + (1/2)T(1I — mm). The tensor Q is diagonal in a 
coordinate system aligned with the principal axes. 



Q 



^-{S + T)/2 ^ 

-{S-T)/2 

y S J 



(5) 



In a frame of reference which is not aligned with the principal axes, Q must be expanded 
in a set of basis tensors T^^ P] which are symmetric, traceless and orthonormal, 

5 

g,;3(x,t) = 5^a,(x,t)T:^. (6) 

i=l 
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Explicitly, these are = Jf'zz , = , /i(x x - y y), = 'x y , = ^2 'x z' 



= y z . 

The Landau-Ginzburg free energy functional F [3] is obtained from a local expansion in 
powers of rotationally invariant combinations of the order parameter Q(x, t), 

Th[Qt] = ^ATrQ2 + l^TrQ^ + hjiTrQ^f + E'{TrCff. (7) 

To this, non-local terms arising from rotationally invariant combinations of gradients of the 
order parameter must be added [3|: 

^e«[<9Q] = ]^Li{da,Ql3^,){do,Qfi^) + ]^L2{daQo,f3){d^Qf3^), (8) 

where a,/3,7 denote the Cartesian directions in the local frame, and Li and L2 are elastic 
constants. Surface terms of the same order in gradients may also be added, but we omit 
them in the present work. 

In the local free energy density Eq.([7]), A = Aq{1 — T/T*), where T* denotes the super- 
cooling transition temperature. From the inequality |(TrQ^)^ > (TrQ^)^, higher powers of 
TrQ^ can be excluded for the description of the uniaxial phase. Thus the uniaxial case is 
described by -E' = whereas E' ^ for the biaxial phase. For nematic rod-like molecules 
B < whereas for disc-like molecules, -B > 0. The quantities C and E' are always taken to 
be positive to ensure stability and boundedness of the free energy in both the isotropic and 
nematic phases. The total free energy including local and gradient terms is 



F = j d^^[^ATrQ'+^BTrQ'+^C{TrQ'f+E\TrQY+hi{d^Q^,)^^^ 

(9) 

The first order isotropic to uniaxial nematic transition at the critical value 5* = 5'c is 
obtained from the equations, 

A = ICS! + ^E'St (10) 
B = -^CS, - 9E'Sl (11) 
The second order uniaxial to biaxial transition at the critical value 5* = 5*^ is obtained from, 

A = , (12) 

45 q 

B = -^E'Sf. (13) 



5 



In the phase diagram of Fig. ([T]), there are hnes separating the nematic phase with S > 
from the discotic phase in addition to the isotropic to nematic transition hnes. There are also 
two uniaxial to biaxial nematic second order lines as well as the two first order isotropic- 
uniaxial nematic lines. One of these transition lines marks the transition to the nematic 
phase with S > and the other for the discotic phase with S < 0. The four phases meet 
at the bicritical Landau point A = B = 0. At the bicritical point, the isotropic to nematic 
transition is second order rather than first order. 

In the absence of thermal fluctuations and hydrodynamic flow, the equation of motion of 
the order parameter can be taken to be relaxational Sj, 

SF 

dtQap{^,t) = -r„^^^— — , (14) 

with 

2 

^al3fii^ = ^[Safj.5i3u + Sai^dp^ — -Saf^S^i^]- (15) 

The kinetic coefficient Tap^v ensures that the dynamics preserves the symmetry and trace- 
lessnes of the order parameter. In the above, F is a constant. 

With the choice of the Landau-de Gennes free energy, the equation of motion takes the 
form 



5t<9a/3(x, t) = -V \{A + CTrQ^)Q^p{^, t) + {B + QE'TrQ') QlJ^, t) 



-L^V^Qapi^, t) - L2 V«(V-,Q^^(x, t)) ]. (16) 

The equations can be projected onto the tensor ial basis T^^ to give a set of equations for 
the independent components ai, 



dta, = -F [{A + CTrQ^)a, + {B + 6E'TrQ^)T^^ Ql^ - L^V^ - T^^^^d^d^a, ]. (17) 

Using the orthonormality of the basis tensors T^/^T^^ = 5ij, this can be expanded to the set 
of five equations. 
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dtai = -r [{A + CTrQ^)ai + {B + 6E'TrQ^)^{al - - + ^ + ^) 

V6 2 2 



6 ' 2 ' ' Vl2 
^5.5.«3)], (18) 

dta2 = -r [(A + CTrQ^)a2 + (S + 6E'TrQ^){-\^a,a2 + ^ - ^) - 

V <j v8 v8 



LiV^as - H^id'y - dl)a, + ^((5^ + d^a^ + 9,5,a4 - dyd.a,))], (19) 
V 12 



9,a3 = -r [{A + CTrQ')as + {B + 6E'TrQ'){-^ + " 

LiV'aa - L2{~^d,dyai + ^{{dl + d^a^ + dyd.a, + d,d,a,))], (20) 

dta, = -r [{A + CTrQ')a, + {B + 6E'TrQ'){^ ^a^^a^^_ 

v6 v2 v2 

LiV^a4 - L2{-^d^dzai + ^{{dl + dDa^ + d:,d^a2 + S^^^^aa + 
V 12 2 



5.5,05))], (21) 
9,a5 = -r U + CTrQ^)a, + (B + ^E'TrQ^)(^ + ^ _ ^) _ 

1 1 
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9.9,04))]. (22) 

Note that these equations are parabohc, non-hnear, contain anisotropic gradient terms and 
are coupled strongly to each other. The efficient computation of the solutions to these 
equations in a variety of physical situations is described in succeeding sections of this paper. 



III. METHOD OF LINES DISCRETIZATION 



The equations for the relaxational dynamics of the orientational order parameter pre- 
sented in the previous section are a set of five coupled non - linear parabolic differential 
equations. No analytical solutions are available for these equations in general and efficient 
and accurate numerical methods must therefore be sought. Below we present such a method, 
based on the strategy of semi-discretisation, whereby an initial value PDE is discretised only 
in the spatial variables to yield a set of coupled ODEs. These ODE's can be solved using 
powerful general purpose solvers. 
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In the literature, this semi-discretisation strategy goes by the name of the 'method of hnes' 
(MOL) ji, ^, since the solutions are obtained along fixed lines in the space - time plane. To 
illustrate this, consider a parabolic initial-value problem in a single spatial variable x, the 
diffusion equation for a scalar ilj{x,t): 

dtiljix,t) = DV'^ip{x,t). (23) 

The MOL discretisation proceeds by restricting the field ip{x,t) to a set of discrete 'col- 
location' points, Xn = nAx, n = 1,2, . . . , N, spaced Ax apart, and then uses a discrete 
approximation for the Laplacian based on these points. The simplest of these approxima- 
tions is based on local polynomial interpolation resulting in finite difference schemes jsl- 
However, the MOL is not restricted to finite differences. When high accuracy is needed, 
global interpolation based on trigonometric or Chebyshev polynomials can be used to gener- 
ate spectral approximations of derivatives [s^. When conservation laws need to be respected, 
finite volume approximations to the derivatives can be used j^. In complicated geometries, 
a finite-element discretisation may be the most appropriate. In this example, the simplest 
approximation is based on nearest-neighbour finite differences, for which 

{V^ij){Xn) = — ^(^(x„+i) - 2tlj{Xn) + H^n-l)) + 0{{Ax)^). (24) 

Inserting this into the diffusion equation, we obtain a set of coupled ordinary differential 
equations 

dttp{Xn, t) = , {^{Xn+l) - 24j{Xn) + ^/'(x„_i)). (25) 

This coupled set of ordinary differential equations, together with initial and boundary con- 
ditions, can now be integrated with a suitable numerical integration scheme which ensures 
accuracy, efficiency, and stability. It is at this stage that the flexibility of the MOL is most 
apparent, since any number of numerical integration schemes can be implemented without 
affecting the spatial discretisation. Depending on the nature of the system of ODEs, the 
optimal choice may be either an explicit scheme, an implicit scheme, or one which is designed 
to handle stiffness. In the above example, for instance, it is well known that an explicit Eu- 
ler integration scheme leads to an instability unless the time step At is constrained by the 

Courant-Friedrichs-Lewy condition At < [AxY /2D 5[. In contrast, an implicit integration 

I 

scheme based on the trapezoidal rule gives the stable Crank-Nicolson update [5|]. In general. 
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semi-discretisation followed by numerical quadrature provides an elegant way of deriving 
many of the well-known schemes for parabolic PDEs. 

The practical implementation of PDE solvers using the MOL discretisation is simple 
since the spatially discretised ODE system can be passed directly to a general purpose ODE 
solver. The complexity, both algorithmic and computational, in the numerical integration 
can thereby be transferred directly to the ODE solver library. 

We have followed the MOL discretisation to construct a solver for PDEs describing the 
dynamics of the orientational order parameter presented in the previous section. We have 
used standard nearest-neighbour second-order accurate finite difference formulae for first and 
second derivatives in the spatial discretisation [10]. More accurate difference approximations 
to derivatives can be obtained using Fornberg's general formula llj. For the temporal 
integration, we have used the ODE solver routines in the GNU Scientific Library. For 
the problems we study here, we find that an explicit multi-stage integrator gives a good 
compromise between accuracy, stability and computation expense. 

Finally, it should be mentioned that the MOL discretisation is not restricted to parabolic 
initial value problems, but also applies to hyperbolic and mixed parabolic - hyperbolic PDEs 
^. This allows the method described in this paper to be extended to situations where effects 
of advection from hydrodynamic flow need to be accounted for. 



IV. APPLICATIONS 



In this section we report benchmarking results using the methodology described in the 
previous section for the following experimentally relevant situations: the properties of the 
isotropic-nematic interface at coexistence, the kinetics of droplet evolution between the 
binodal and spinodal lines and spinodal decomposition into the nematic phase. In each 
case, we compare our numerical results with available analytical and numerical results in 
the literature. 

A suitable non-dimensionalisation of the equations of motion before they are discre- 
tised is essential for controlling artefacts and error introduced by the discretisation. For 
our problem, we non-dimensionalise the governing equations ( IT6l) using the values of di- 
mensionful quantities at coexistence (T = Tc). The resulting dimensionless quantities are 
superscripted with an asterisk, while the dimensionful quantities are subscripted with c. 
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We non-dimensionalise the order parameter Q*^ = Qa/s/Sc by the value of the strength 
of ordering at coexistence Sc = — 25/9C. The non-dimensionahsed strength of ordering is 
then S* = S/Sc = + a/1 — 2AAC/B'^)/A. The gradient terms in the free energy are non- 
dimensionahsed by the length Ic = a/54C(Li + 2L2/3)/i?^, which is related to the width 
of the isotropic-nematic interface at coexistence. The free energy is non-dimensionalised as 
F* = F/Fc, where F^. = 9CS^/16 is proportional to the free energy barrier between the 
isotropic and nematic minima at coexistence. Finally, time is non-dimensionalised by the 
characteristic relaxation time r = (TFc)^^. 

To solve the equations of motion Eq.( fT7|) . we choose the discretisation length Ax = 
Ay = 1 and the integration time-step At = 1. These define simulation units of length 
and time. To ensure that discretisation errors and artefacts are kept to a minimum, the 
discretisation length must be much smaller than the characteristic length Ax = Ay <^ Z^- 
The integration time step must be much smaller than the characteristic time scale At ^ r. 
The dimensionless discretisation scales must then satisfy Ax* = Ax/lc ^ 1, Ay* = Ay/lc -C 
1 and At* = At/r ^ 1. Our simulations are performed maintaining the above conditions, 
on grid sizes ranging from 32 x 32 to 256 x 256, with periodic boundary conditions, using 
the 4th-order Runge - Kutta method as implemented in the GNU Scientific Library for the 
ODE solver. 



A. The isotropic-nematic interface at coexistence 

The isotropic-nematic interface was studied in a seminal paper by de Gennes within the 
framework of a Landau- Ginzburg description. To render the problem analytically tractable, 
de Gennes made a specific assumption regarding the variation of the components of the order 
parameter across the interface. For an infinitely extended interface where, by homogeneity, 
variations perpendicular to the interface alone are allowed, de Gennes assumed that the 
only quantity which changed across the interface was the strength of ordering 5*. In the de 
Gennes ansatz, there is no biaxiality and no variation of the director across the interface. 
This reduces a problem with five degrees of freedom to a more manageable problem involving 
only a single degree of freedom. The variation of the ordering strength 5* along the coordinate 
z normal to the interface located at zq can then be obtained analytically as, 

^(^) = ^(1-tanh^^), (26) 
10 



where w = \/2/ S* is the non-dimensional interfacial width. 

We have verified this remarkable ansatz with a direct numerical solution. In our numerical 
calculations, a strip of nematic interface was sandwiched between two isotropic domains with 
periodic boundary conditions. The system was then allowed to relax to the minimum of the 
free energy. The parameters were chosen such that the width w ^ Az = 1, ensuring that 
discretisation errors were kept to a minimum. 

The resulting profiles for the variation of S and T are shown in Fig. ([2]). The values 
obtained for T are consistent with de Gennes' assumption of vanishing biaxiality. The 
variation of S at each of the two isotropic-nematic interfaces were fitted, using the least- 
squares method, to the analytical profile, with saturation value of the order Sc, the location 
of the interface zo and the interface width w as fitting parameters. As shown in the inset 
to Fig. ([2]), fitted values of w agree remarkably well with the analytic result for a range 
of parameter values. The agreement is accurate to within a fraction of a percent. Similar 
results were obtained for the saturation value of the order parameter. This benchmark 
clearly demonstrates the accuracy of the MOL scheme in reproducing the equilibrium limit 
of Eqn.dll). 

The de Gennes ansatz also predicts that the energies of planar and homeotropic anchoring 
are identical when elasticity is isotropic {L2 = 0). We compared the value of the free energies 
of the interface for both planar and homeotropic anchoring of the director. To machine 
precision, these answers are identical. Our results for this problem represent the first direct 
verification of the de Gennes ansatz retaining all degrees of freedom of the orientational 
tensor. 

B. Nematic droplet in an isotropic background 

Since the isotropic-nematic transition is first order, there exists a regime of parameters 
where the kinetics proceeds by nucleation. In the phase diagram of Fig. ([T]) this regime 
is bounded by the binodal and spinodal lines. A droplet of nematic phase in an isotropic 
background shrinks if it is smaller than a critical size Re- Droplets which are larger than 
Rc grow in size till they expand to the size of the system. A droplet of size Rc is thus 
metastable. 

The development in time of a fiuctuation-induced droplet is a key factor in determining 
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the nucleation rate [6j. Unlike conventional isotropic fluids, the nucleus in a nematic is not 
expected to be spherical 1^. Allowing the shape to deviate from perfect sphericity reduces 
the total energy, which is a sum of the elastic energy associated with the director deformation 
in the bulk and a surface energy associated with the anchoring condition at the interface of 
the droplet. Thus, determining the droplet shape of least energy involves a minimisation 
over the strength of ordering as well as the director degrees of freedom 

In our simulations we have studied the time evolution of droplet shapes. For studying 
conformational changes, our initial condition was chosen to be S* = Sc, with the director 
anchored at an angle of 7r/4 with the nematic - isotropic interface inside the droplet. We 
took S* = in the isotropic region. We relaxed the system at a temperature intermediate 
between the binodal and spinodal temperatures. 

For droplets larger than the critical radius, the nematic region was observed to grow in 
size, flnally evolving into the full nematic state. For the parameters used, the alignment of 
the director did not change significantly during the evolution of the droplet. Figs. ([3]) shows 
three stages in the evolution of a single, initially circular droplet, illustrating how the droplet 
shape evolves into an elongated ellipsoidal form. In the absence of elastic anisotropy L2 = 0, 
the droplet remains circular. In the presence of elastic anisotropy, the droplet orients along 
the direction of nematic order for planar anchoring (L2 > 0) and perpendicular to it for 
homeotropic anchoring (L2 < 0). 

In our numerics, we quantify this geometrical change through measurements of the change 
in the aspect ratio. The aspect ratio, defined as the ratio of the major to minor axis of the 
ellipse, are indicated in Figs. ([3]). This is calculated by extracting a contour at a fixed value 
of S (say Sc/2) and fitting it with an ellipse. 

To the best of our knowledge, there are no analytical expressions for order parameter 
variations for nematic droplets. Thus, our results cannot be compared directly against 
analytical theory. However, the results obtained are in qualitative agreement with previous 
work on the shape of nematic droplets as a function of anchoring strength jisl. 15 . 16|. 



C. Spinodal coarsening 

In the previous section, we studied the dynamics of droplets when the quench was to 
a temperature between the binodal and the spinodal. For a quench below the spinodal 
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temperature, the isotropic phase becomes locally unstable to a nematic perturbation and 
the system proceeds to the nematic phase by spinodal decomposition. Coherent regions of 
local nematic order develop in time, with a distinct axis of order in each of these domains. 
Topological defects form at the intersections of these differently ordered domains. Coarsening 
proceeds through the annihilation of topological defects, increasing the correlation length of 
local orientational order. 

To study the coarsening kinetics, we start from a random initial configuration where 
we draw the order parameters S and T randomly from a Gaussian random variable with 
variance proportional to Fc, ensuring that < T < S. We obtain cos 6' from an uniform 
distribution between -1 and 1, and choose similarly between and 2tt to generate the 
director, codirector and the joint normal. We then relax the system from this initial condition 
at a temperature below the supercooling spinodal temperature. The data presented below is 
averaged over 100 different initial conditions for a 256 x 256 system with periodic boundary 
conditions. 

From the coarsening simulations we obtain the strength of ordering, the biaxiality and 
the director. The director is used to construct the schlieren plots shown in Figs. (HI). These 
plots are constructed by first projecting the director into the x — y plane, finding the angle x 
made by this projection with an arbitrary axis (say the x- axis) and then computing sin^(2x). 
The presence of both integer and half-integer defects is clearly visible in these plots as the 
meeting points of four and two dark brushes, respectively. In the corresponding plots for 
the strength of ordering, the defects are clearly visible as localised regions where S rapidly 
decreases. This is the core region of the topological defect, shown in Figs. ([5]). We confirm 
the surprising finding that there is strong biaxial ordering inside the defect core [l3] • These 
results are in perfect qualitative agreement with both theoretical predictions and previous 
numerical results 17l |. 



To make a quantitative comparison with previous work, we compare results for the time- 
development of correlation functions during coarsening. We calculate the real-space corre- 
lation function, 



and its Fourier transform, 



"^''''^ /rf3kg.^(k,t)g^«(-k,t)- ^'''> 
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Theoretical predictions and analytical work have verified that these correlation functions 



have a scaling form C{r,t) = J^[r/L{t)], and S{k,t) = L'^{t)g[kL{t)] [1$. Here, C{r,t) = 
J2\r\=r^(^'^) S{k,t) = X]|k|=fc '^C^' ^) angular averages of the correlation functions 
in real and Fourier space respectively. The length L{t) is extracted from the real-space 
correlation function using the implicit condition C{r = L(t),t) = 1/2. In Fig. (jS]) we confirm 
that the real-space correlation function does indeed scale as expected. Our numerical data 
for the scaling function is in close agreement with an analytical calculation for the n = 2 
model due to Bray and Puri |19]. In Fig. ([7]) we show the corresponding scaling of the 
Fourier space correlation function. The wavenumber {k) is the root of the second moment 
of the S(k, t) defined by 

{k)' = ^^ = 11 k'S{k, t)l J2 S{K t). (29) 

The inset shows the growth of the length scale as a function of time. Theoretically, this is 
expected to grow as a power L{t) ~ t". Our estimate for this exponent is a = 0.5 ± 0.005. 

Our results are consistent with both analytical predictions and an earlier numerical sim- 
ulation. The Fourier space correlation function is expected to exhibit a short-wavelength 
scaling S{k,t) ~ k~^ known as a generalised Porod's law |l^. We see a clear range of 
wavenumbers where the Porod scaling is obtained. At very short wavelengths, correspond- 
ing to the size of the defect core, the Porod scaling breaks down. We see evidence for this 
as well, where the very highest wavenumbers in Fig. ([7]) show deviations from the Porod 
scaling. Our numerical results for spinodal decomposition, then, agree both qualitatively 
and quantitatively with theoretical results and previous numerical work [20^. There are 



numerical results reported in the literature are confiicting 2l|, |22|. We indicate possible 
reasons for this discrepancy in the following section. 

D. Discussion 

For the nematic-isotropic interface, elastic anisotropy (L2 7^ 0) induces biaxiality near the 
interface. This problem has been studied using various approximations by several authors. 
Sen and Sullivan [23| introduced a symmetry-adapted parametrisation of the order parameter 



which reduces the independent degrees of freedom to three. Using this parametrisation Popa- 



Nita and Sluckin 



24 



25| obtained numerical solutions for the variation of S and T across 
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the interface. The parametrisation used in these calculations is only applicable to the cases 
of planar or homeotropic anchoring at the interface. For a general anchoring condition this 
parametrisation breaks down, as is obvious from the fact that the order parameter has, in 
general, five independent components. 

With the present method, the same problem can be studied without any approximation, 
retaining all degrees of freedom of the orientation tensor. This allows us to obtain the 
variation of the order parameter in situations where the symmetry conditions of Sen and 
Sullivan do not apply. These include curved interfaces and the effect of higher order elasticity. 

With regard to the problem of nematic droplets a detailed study becomes possible, espe- 
cially in three dimensions. Indeed the only numerical work we know of is the Monte Carlo 
simulation of Cuetos and Dijkstra [6]. The advantage of the Landau-de Gennes approach 
for this problem is that the free energy of anchoring need not be postulated, (as is done 
when only the director degrees of freedom are retained and a Rapini-Papoular 26|] free en- 
ergy used), but is effectively present in the Landau-de Gennes free energy itself. This is 
the calculational strategy used by Prinsen and van der Schoot [l^. Non-trivial director 
configurations have been proposed for different strengths of bulk-to-surface coupling. Such 
predictions can be verified cleanly with the present method. 

For nematic coarsening in three dimensions, simulations within the Landau-de Gennes 
framework have not been performed so far. The extended nature of topological defects, 
which are lines rather than points in three-dimensional nematics, makes the dynamics of 
the coarsening problem quite different from two dimensions. The strength of the Landau- 
de Gennes approach is particularly clear here. Having access to the strength of ordering 
makes it easy to locate the defect lines as tubes of zeros of the ordering strength. This is 
considerably easier than using a Burgers-like circuit integral of the director field to locate 
defects as is done in the work of Zapotocky et al j21 |. 

The dynamics of defect lines has implications in cosmology, where the Kibble mechanism 
271] predicts a scaling relation for the number density of defects. Experiments studying the 
Kibble mechanism in liquid crystals have been performed 28|]. However, we know of no 
numerical study of the Kibble mechanism along the lines proposed here for nematic liquid 
crystals. This provides further impetus for three-dimensional simulations. 

There are two alternatives to the MOL discretisation for tensor order parameter descrip- 
tions of the nematic phase. These are the cell dynamical scheme of Oono and Puri 291] as 
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22 1, and the lattice Boltzmann 



implemented by Zapotocky et al 2l|] and Dutta and Roy 

method of Denniston et al 30^. Our method differs in an important way from both these 
approaches, in that it provides a direct discretisation of the governing equations of motions. 
In cell-dynamical simulations, a local map is constructed whose fixed points coincide with 
the extrema of the local part of the free energy. This local map is phenomenological and 
does not follow from a Landau expansion. The advantage of this method is that it produces 
sharp interfaces, ensuring a rapid scale separation between the microscopic interfacial length 
and the macroscopic domain size. In coarsening simulations, this reduces the non-universal 
offset time at which the system enters the scaling regime. 

In comparison, the Landau-de Gennes free energy leads to a less sharp interface. Conse- 
quently, the structure of the core region of topological defects is different in the two methods, 
being smaller in cell dynamics and larger in the present method. We see evidence of this 
in the dynamic structure factor. Fig. ([7]), where there is a violation of Porod scaling at the 
highest wavenumbers. As shown in Figs. ([5]) the core region of our defect spans several lattice 
spacings in each direction. Thus, as regards scale separation between what are microscopic 
and macroscopic lengths, the cell dynamical scheme fares better. 

On the other hand, when physics at the scale of the interface is important, as in the 
problem of the planar interface and nematic droplets, the MOL discretisation offers a clear 
advantage. Due to the special nature of the cell dynamical Laplacian, it is not clear how 
to generalise the cell dynamical method to include elastic anisotropy (corresponding to the 
gradient term with coefficient L2) or higher order terms involving antisymmetric contrac- 
tions. Due to the phenomenological nature of the cell dynamical map, it is difficult to make 
direct contact with experiment, a procedure which is straightforward in the MOL when the 
non-dimensionalisation we outlined above is applied. 

In the lattice Boltzmann method, the governing parabolic equations are replaced by a 
hyperbolic superset 3]| through the introduction of a distribution function for the tensor 
order parameter. As with cell dynamics, the Landau-de Gennes equations are not solved 
directly. Rather, lattice Boltzmann relies on a temporal scale separation which allows the 
hyperbolic equations to mimic parabolic behaviour. The lattice Boltzmann method needs 



to be carefully formulated [32 
appears to be non-universal 



_to avoid the presence of micro-currents and this procedure 
We find no evidence of spurious micro-currents in the 



MOL discretisation presented here. The lattice Boltzmann method includes the coupling of 
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order parameter and flow. It is entirely straightforward to extend the MOL discretisation to 
include order parameter advection. To include coupling to the fluid momentum, we advocate 
a hybrid strategy, where the order parameter dynamics including advection is solved by a 
MOL discretisation, but the dynamics of the fluid momentum including order parameter 
stresses is solved by the lattice Boltzmann method. 

In terms of algorithmic complexity, cell dynamics, lattice Boltzmann and the MOL with 
flnite-differences appear to be matched, since both use fairly local information to calculate 
derivatives and involve explicit temporal updates. For N degres of freedom the algorithmic 
complexity of all three algorithms is 0{N). On the other hand, the storage requirements 
of the lattice Boltzmann formulation are larger by a factor of 6 to 9 in two-dimensions, 
and about 15 to 19 in three dimensions. On an Intel(R) Core(TM)2 Duo CPU with speed 
2.66GHz, our code takes 0.33 seconds for one time step once a lattice size of 256 x 256. 

V. CONCLUSION 

This paper has proposed an efficient numerical scheme based on the method of lines for the 
solving the Landau-de Gennes equations of nematodynamics. The numerical results obtained 

in previous sections are in excellent agreement with analytical results where available and 
consistent with previous numerical data. We expect that the method presented here will 
find broad application in exploring the rich physics of the nematic phase of liquid crystals. 
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FIG. 1: (Color online) The mean field phase diagram obtained from the Landau-de Gennes free 
energy ([7]). The phase boundaries are given by solutions of the following algebraic equations : 
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the supercoohng spinodal; gB'^E' + 8B'^{C^ - 36ACE') = 192A{C'^ - AAE' f for the superheating 
spinodal. 
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FIG. 2: (Color online) Variation of the degree of alignment (S) and biaxiality (T) across the 
nematic - isotropic interface with planar anchoring. Symbols represent numerical data, solid curves 
are the de Gennes ansatz ()26p . The numerical parameters are A = 0.0035, B = -0.5, C = 2.67, 
E' = 0,Li = 0.01,^2 = 0,r = 1/20, grid size 8 x 512 and nematic strip width = 256. The inset 
shows the variation of the interfacial width with the elastic constant Li. The expected quadratic 
variation is accurately reproduced. 
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FIG. 3: (Color online) Time evolution of the degree of alignment in the relaxation of a nematic 
droplet. The initially circular droplet develops an anisotropy and becomes elliptical, the two- 



dimensional analogue of the ellipsoidal droplets (factoids) seen in three dimensions. Figs. 3(a) 

3(e) and 3(i) shows the circular droplet at time t = 0. The time evolution with L2 = of the 

circular droplet are shown in Figs. |3(b)||3(c) and |3(d)] The elliptical conformations with L2 > are 

shown in Figs. |3(f) 3(g) and |3(h)| These have aspect ratio 1.1191, 1.3046 and 1.5037 respectively. 

With L2 < conformations are shown in Figs. 3(j), 3(k)| and |3(l)[ These have aspect ratio 1.0122, 

1.0718 and 1.1316 respectively. The numerical parameters are A = 0.001, B = -0.5, C = 2.67, 

E' = 0, Li = 0.0236, r = 1.0, grid size 128 x 128 and the droplet radius = 20. The elastic constant 

L2 = lOLi for the L2 > and L2 = -Li for L2 < 0. 
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(d) (e) (f) 

FIG. 4: (Color online) Schlieren textures in a coarsening nematic. Topological defects of both 
integer charge (top right in (c) and (d)) and half - integer charge (throughout) are observed in the 
simulations, (a) shows the random configuration at time t = 0. (b), (c), (d), (e) and (f) show the 
dynamics at different times t = 10'^, 5 x 10'^, 2 x 10^,6 x 10^ and 10^ respectively. The parameters 
chosen are A = -0.1, B = -0.5, C = 2.67, E' = 0.0, Li = 1.0, L2 = and T = 1.0, grid size 256 x 256 
for 10^ time steps. 
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FIG. 5: (Color online) The schlieren texture around a half integer charge defect is shown in (a), 
(b) shows the density plot of the uniaxial order parameter, while (c) exhibits the variation of 
the uniaxial and biaxial order parameter alongj^ line passing through the defect core. Note the 
presence of strong biaxiality within the defect core as seen previously The sharpness of the 




FIG. 6: (Color online) Data collapse of the direct correlation function C(r) with scaled distance 



Mi for the 0(2) model 



r/L{t) for different times. The symbol A depicts the Bray - Puri function 
: fBp{x) = B'^{0.5,1.5)F[0.5,0.5,2;exp{—x'^)]exp{—x^/2)/TT. The inset shows the unsealed corre- 
lation function at different times. Numerical parameters chosen were A = —0.1,B = — 0.5;C = 
2.67; E' = 0, Li = 1.0, L2 = 0, T = 1/20 on a 256 x 256 grid 
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FIG. 7: (Color online) Data collapse of the structure function S{k,t) at different times. The dash- 
dot line has a slope of —4 indicating the validity of generalised Porod's law for 0(n) systems. The 
departure from Porod's law at high k is due to the finite core size of the defects as discussed in the 
text. The inset shows the time dependence of the correlation length L(t). The length grows as a 
power law with an exponent of 0.5. The maximum value of the correlation length is approximately 
1/4-th the system size, ensuring the absence of finite-size artefacts. Numerical parameters chosen 
are A = -0.025, B = -0.5, C = 2.67, E' = 0, Li = 0.1, La = 0, F = 1/20 on a 256 x 256 grid. 
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